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Piecewise smooth maps are known to exhibit a wide range of dynamical features including nu¬ 
merous types of periodic orbits. Predicting regions in parameter space where such periodic orbits 
might occur and determining their stability is crucial to characterize the dynamics of the system. 
However, obtaining the conditions of existence and stability of these periodic orbits generally use 
brute force methods which require successive application of the iterative map on a starting point. 

In this article, we propose a faster and more elegant way of obtaining those conditions without 
iterating the complete map. The method revolves around direct computation of higher powers of 
matrices without computing the lower ones and is applicable on any dimension of the phase space. 

In the later part of the article, we compare the speed of the proposed method with the other popu¬ 
lar algorithms which shows the effectiveness of the proposed method in higher dimensions. We also 
illustrate the use of this method in computing the regions of existence and stability of a particular 
class of periodic orbits in three dimensions. 
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I. INTRODUCTION 

Presence of stable and unstable periodic orbits play 
a vital role in determining the dynamical properties of 
a system. While presence of stable periodic orbits form 
the basins of attraction, the unstable periodic orbits play 
a crucial role in forming chaotic orbits and the basin 
boundaries of the attractors. In fact appearance, dis¬ 
appearance or change in stability of these periodic orbits 
are responsible for most bifurcation phenomena. Hence 
an important component in characterising a given sys¬ 
tem is determining the parameters for which a particular 
periodic orbit might exist. 

Piecewise smooth maps are useful in describing sys¬ 
tems whose evolution is given by different smooth func¬ 
tions in different regions of space. Examples of piece- 
wise smooth systems include switching electrical circuits 
[1, 2], impacting mechanical systems [3], walking robots 
[4, 5], cardiac dynamics [6] and neural spiking [7]. Apart 
from their applicability, piecewise smooth maps have also 
attracted attention due to their rich dynamical proper¬ 
ties. Even maps which are piecewise linear and have 
non-linearities only across a switching manifold, exhibit 
features like robust chaos [8] and existence of infinitely 
many co-existing attractors [9]. Study of piecewise linear 
maps attains even more importance as they can describe 
the behaviour of the piecewise smooth systems near the 
border separating the partitions of the phase space, in 
which case the matrices take a specific form. This article 
aims to study the conditions of existence and stability 
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of periodic orbits in such piecewise linear maps in any 
arbitrary dimension. 

Periodic orbits in piecewise linear one, two and three 
dimensional maps have been studied [10-12]. Some of the 
works have also computed the regions in parameter space 
where specific types of orbits exist and are stable [12-15]. 
However most of the studies used repeated application of 
the iterative map to obtain them. An iterative technique, 
invented by the Russian mathematician Leonov [15-17] 
has also been used to obtain the existence condition of 
periodic orbits [13, 18, 19]. Moreover, the methods used 
there were specific to the type of orbit being analysed. 

In this work, we seek to obtain a generic and faster 
method of finding the existence and stability criteria for 
a general class of periodic orbits. The technique devel¬ 
oped might be extended to any periodic obit in piecewise 
linear maps of any dimension. Moreover, as the tech¬ 
nique is essentially an algorithm to obtain higher iterates 
of the map without going through the intermediate ones, 
it can be used to simplify extraction of other relevant 
information about the map. Additionally, the maps we 
work on are in their normal forms. Hence the ambit of 
the results obtained spreads across all piecewise smooth 
maps as appropriate coordinate transformation near the 
border of non-smoothness yields the normal form of the 
maps. 

The paper is organized as follows. After defining 
the notations and conventions used in the article in the 
Section-2, we derive the generic conditions for existence 
and stability of periodic orbits in a dimension indepen¬ 
dent form. We then develop a technique which exploits 
the form of the map to obtain the conditions of existence 
and stability in N dimensions. We showcase the utility 
of the technique developed by computing the parameter 
regions in which certain classes of stable periodic orbits 
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exist in a three dimensional piecewise linear map in the 
next section. Thereafter, we take the specific case of two 
dimensional systems, where the simplicity of the map al¬ 
lows us to obtain stronger analytical results. Finally, we 
devote a section to exhibit the computational efficiency 
of the technique developed by comparing it with other 
traditionally used methods for computing the existence 
and stability criteria of periodic orbits in piecewise linear 
maps. 


II. THE PIECEWISE LINEAR MAP IN THE 
NORMAL FORM 

In earlier literature it has been shown that, on proper 
choice of axes, any N dimensional piecewise smooth con¬ 
tinuous map can be linearised near the border to have 
the following form [10] 


and 


GM) = 


where Ml, Mr £ 
form given by [20] 


MlX + ( : X\ < 0 
MrX + £ : X\ > 0 


( 1 ) 


rxN are matrices in their normal 
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with d\ J ^ being the coefficient of A* in the characteristic 
polynomial of Mj\ X = (xi,..., Xn) T £ being a 
generic point in the phase space and ([ = (p, ...,0) T . 
Since the coefficient of A* in the characteristic polynomial 
of a matrix is the sum of eigenvalues of the matrix taken 
j at a time (apart from an alternating sign); we recast 
(2) as 




(J) 


-d 

\-d 


(j) 
N -1 
(J) 

N 


1 0 
0 1 

0 0 
0 0 


o\ 

0 

1 

0 ) 


: J £ {L, R} (2) 


Mj = 


/ ( J) 

( Pi 

( J ) 

-P\ 


(—l )^ 2 P$U 

V (~l ) N ~ 2 pW 


1 0 
0 1 

0 0 
0 0 

AJ) 


o\ 

0 

1 

0 / 


Je{L,R} 


( 3 ) 

for further analysis. Here p\ J 1 is the sum of the eigenval¬ 
ues of Mj taken i at a time, i.e., p\ JS> = (—1 ) l ~ 1 d[ J \ For 
instance, in two and three dimensions, the matrix in (3) 
takes the form, 


Mj = 


Tj 1 
-Sj 0 


: J £ {L, R} 


( 4 ) 



Je{L,R} 


( 5 ) 


respectively. Here tj and Sj are the trace and determi¬ 
nant the Mj; and aj is the sum of the eigenvalues of Mj 
taken two at a time. 

Any p-periodic orbit in such a system can be repre¬ 
sented by a sequence of points {Xo,..., Xi ,..., A p _i} 
such that Xi + i = G^Xi) where i is a natural number or 
zero, and X p = Xo. We can also associate a symbol to 
each point on the periodic orbit depending on the par¬ 
tition in which it lies. If a point of the periodic orbit 
has Xi < 0, it is assigned the symbol L (meaning left). 
Otherwise it is assigned the symbol R (meaning right). 
This associates a sequence of R and L to each periodic 
orbit. For example, an LLLR (also written as L 3 R) orbit 
consists of three points with Xi < 0 and a single point 
with Xi > 0. To remove the ambiguity arising from pos¬ 
sible cyclic permutation of the points, in this article we 
adopt the following convention: We number the points 
Xi such that Xo is assigned the symbol R and X p is 
assigned the symbol L. However while referring of the 
orbits, we collect the symbols in the reverse order such 
that the symbol for the orbit starts with L and ends with 
R. In this article, for sake of simplicity, we would restrict 
our explicit analysis to orbits of the form L m R n where 
to and n are natural numbers. However, the theory de¬ 
veloped here can be extended to the general finite period 
orbit L ni R n2 .. .L ns ~ 1 R nB . Ways to extend the theory 
for the general cases will be indicated wherever required. 
Note that since the map in each of the partitions is lin¬ 
ear, a periodic orbit lying completely within any one of 
the partitions does not exist. 


III. EXISTENCE AND STABILITY OF 
PERIODIC ORBITS 


In 1959, Leonov [16, 17] did a detailed study of 
nested period adding bifurcation structure occurring in 
piecewise-linear discontinuous ID maps. The algorith¬ 
mic way proposed in his work was recently used[13] to 
analyse border collision bifurcations in one dimensional 
maps. We extend the analysis to obtain the existence 
criteria for period orbits in N dimensions. 

Consider an L m R n orbit formed by the points 
{Xo, ■ ■ ■, X m+n _i}. According to the conven¬ 
tion described in the last section, we assume the 
points Xo,-..,X m _i have Xi > 0 and the points 
X m ,..., X m+n _i have Xi < 0. Hence the evolution of 
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Xq under the map G M is given as, 

X! = M r X o + C 

X 2 = M%X 0 + (/ + Mr) C 

X 3 = MrXq + (/ + Mr + Mr) £ 


is to compute powers of Ml and Mr. Typically this 
is done algorithmically by brute force matrix multiplica¬ 
tion. The following sections aims to provide an algebraic 
way to compute the powers of the matrices and hence to 
analytically compute the regions of existence and stabil¬ 
ity of the periodic orbits. 


X n -i =M£- 1 Xo + <^, n _iC 
X n = M%X 0 + c 

X n+ i = M l MrXq + M L cf> R , n C + c 


IV. COMPUTING M n FOR N DIMENSIONS 

In this section, we present a technique for computing 
n th powers of the matrix M , given by 


X m+n -i = 0 + (M™~ l <j> R , n + C 

X m+n = M™ MrXq + + (t>L,m) C- 

Here / is the N x N identity matrix and <f>j t k = 1 + 
Mj + M 2 j+ ... + Mj- 1 . By using the formula for GP of 
matrices, it can be written as, 


f Pi 10 
P‘2 0 1 


M = 


(-1 )*- 2 PN-1 0 0 

V PN 0 0 


°\ 
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( 10 ) 


= (I - Mj)- 1 (/ — Mj) :JG {L, R} (6) 
if I — Mj is invertible. 

On substituting X m+n = Xq, we get the expression for 
Xq in terms of known quantities, 


Note that this matrix is the same as the matrices Ml 
or Mr defined in (3) except that the subscripts L and R 
are ignored. This is done for notational simplicity. The 
results derived here are applicable to both Ml and Mr. 

The n th power of any general N x N matrix 


Xq = (I - M^Mr)- 1 (M?cj>R in + ct>L,m) C (7) 

when (/ — is invertible. Once Xq is determined, 

all other points of the orbit can be calculated by evolving 
Xq under the map G M . For existence of the L m R n orbit, 
we need to ensure that all points of the periodic orbit are 
in their correct partitions. 

For periodic orbits to be stable, we require the trace T 
and determinant A of the ordered product of the Jaco¬ 
bian matrices at each point of the periodic orbit to follow 

(1- A) <T< (1 + A). (8) 
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can be written in the form 

/V") .. 
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n(n) 
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( 11 ) 


( 12 ) 


Since the Jacobian at any point with x\ < 0 is Ml and 
that at any point with x\ > 0 is Mr, the matrix whose 
trace and determinant are in question is 

For the generic case of L ni R n2 ... L ns - 1 R n, ‘ orbits, the 
algorithm of finding the existence and stability conditions 
of the orbits remain the same. On explicit calculation, 
the expression for the the point Xq on the periodic orbit 
comes out to be 


X 0 = 


-1 


7 -n m k 


1 


e n m "k 



(9) 

where K is the symbol L or R when j is odd or even 
respectively. The condition for stability remains identical 
to (8), with the only difference that now T and A denote 


the trace and determinant of the matrix n 

3= 1 

A close look at equations (7), (8) and (9) reveal that, 
the primary requirement in evaluating those expressions 


The sequence of matrices A, A 2 , A 3 ■ ■ ■ A n is composed of 
sequences of its individual elements, and therefore N 2 
such independent sequences are needed to construct the 
matrix A n . However, we show that due to the structure 
of the matrix M in (10), only N sequences are required 
to construct the matrix M n . 

Let each of the N sequences be denoted by Ti. Let 
Tij denote the j th term of the sequence Tj. Also note 
that for notational simplicity in the further analysis, we 
start numbering the terms of the sequence from 2 — N 
instead of 1. For example, in order to compute M n 
in three dimensions, we would require three sequences: 
Ti, r 2 and Ts; and the terms of the sequences will be 
numbered as r^!, Fi 0 , T^i, r 12 , • • • for the first se¬ 
quence, r 2j _i, r 2> o, r 2j i, r 2j2 , • • • for the second sequence 
and r 3 o, Tsp, T 3 i2 , • • • for the third sequence. 

The terms of this sequence shall be defined iteratively. 
However, before giving the iterative relation, we first give 
the first N terms of the sequences. This initialisation is 
done such that the N elements of the i th row of M give 
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the first N terms of in the reverse order, i.e., 


/ r bl 

Ti.o • 

• Ti , 2 _JV^ 

r 2 ,i 

r 2 , 0 • 

• r 2j2 _Ar 

V^AT,! 

r./v,o • 

• Fn^-nJ 


Again taking the three dimensional case as an example, 
we would have the initialisation as, 

Ti.-i = 0, ri l0 = IjTip = r 

1 2.-1 = i) r 2i0 = o,r 2 ,i = —a (14) 
r 3 ,-i = o, r 3>0 = o,r 3 ,i = s. 

Note the due to the initialisation given above, the 
terms from rj i2 _Ar to r^i are defined. Now, the further 
terms of the sequence are defined iteratively as 

N 

Tij = P k : 3 G [2, oo). (15) 


1.0 



T,. 


With this definition of Fjj, the matrix M n can be 
written as 


^i,n I\ n _i ... r 1 )Tl _(jv_i)\ 


M n = 


2,n-(N-l) 


(16) 


yT7V,n bjv.n-l N,n—(N—1) J 

or more explicitly as 


[M"] • ■ = Vn e N. (17) 

The proof of the result hinges on the structure of the 
matrix M. Note that (17) is bound to be satisfied for n = 
1 due to the way the series Fj is initialised in (13). Also 
note that, apart from the first column, the only non-zero 
elements of the matrix M are in the superdiagonal and 
are all equal to 1. Hence, while obtaining M k+l from M k , 
when M is multiplied from the right, the columns of M k 
are simply shifted to right, except the rightmost column 
which is lost. Therefore, only the first column needs to 
be computed, which is given by (15). The detailed proof 
of the result is given in Appendix A. 

Hence in order to compute M k+1 from M k , one needs 
to compute only N new elements corresponding to the 
first column of the matrix; as compared to computing 
N 2 new elements for a generic N x N matrix. 

Having described the algorithm of finding the n th 
power of an N dimensional matrix, we apply the algo¬ 
rithm to 3 dimensional piecewise linear maps to compute 
the regions in parameter space where stable periodic or¬ 
bits might exist. 


V. APPLICATION: STABLE L n R ORBITS IN 3 
DIMENSIONS 

In 2012, parameter regions where stable L n R orbits 
exist in two dimensions were computed to find the pa¬ 
rameter regions of multiple attractor bifurcations [12]. 


FIG. 1: Regions of the projected parameter space where 
stable L n R orbits exist. The plot has been made by 
setting o\l = ctr = 1.4, 5l = Sr = 0.7 and n = 1. 


In this section, we extend the result to demonstrate the 
use of the technique developed in this article to find the 
regions of existence of stable L n R periodic orbits in 3 
dimensions. As the parameter space is 6 dimensional, 
we show the two dimensional projection of the plausible 
regions in parameter space. 

In three dimensions, the normal matrix (10) takes the 
form 

/r ! 0\ 

M = —a Oil. (18) 

\ 5 0 0 / 

Hence, by (16) 

ff'l.n I\ n _i Ti iTI _ 2 \ 

M n = r 2 , n r 2 , n -i r 2 ,„_ 2 ( 19 ) 

\r 3 ,n r 3n _i r 3 „_ 2 y 

where 

r i,n = rr- oT i)n _2 + (5r ijn _ 3 (20) 

for * € {1, 2, 3} and n > 1. For —1 < n < 1, the terms are 
taken from M according to (14). Using these expressions, 
we compute the required powers of Mr and Mr, and 
substitute them in (7) and (8) to obtain the required 
conditions for existence and stability of L n R orbits for 
various values of n. The parameter values for which the 
stable orbits exist are shown in Fig. 1 where the plausible 
regions are shown in the two dimensional projection of 
the six dimensional space. 

Now let us consider the special case of 2 dimensions 
where further simplification can be done and the n th 
power of M can be computed non-iteratively. 
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VI. COMPUTING M" FOR 2 DIMENSIONS 

In two dimensions the matrix appearing in the normal 
form map is given by (4). Dropping the subscripts for 
notational simplicity as in the previous section gives us 


which can be recast into the form of (24) with the defi¬ 
nition of a n as 


a„ 


A”"*" 1 — A” +1 

Ai — A2 


(28) 


m =(_ t { ;) ■ (2D 

Prom the results obtained in the previous section, it 
can be said that M n would be determined by two inde¬ 
pendent sequences a and b and would be of the form 


M "=(£ tS) ■ (22 > 

However on explicit calculation (as done in Appendix 
B) it can be shown that the terms of the two sequences 
are related as 


bi = —Sen -1 Vi > 0. (23) 


The detailed proof of this result is given in the Appendix 
D. It can also be directly seen that (28) satisfies (25). 

Hence in order to compute M n for a 2 x 2 matrix in 
the form (21), we use (26) or (28) to get a n , a n -i, a n -2 
and form the matrix given in (24). 

The form of the matrix M n in (24) may be substituted 
directly into (6) to obtain the sum of the GP as 


(fin 


( fn fn — 1 

\-Sfn -1 1 - Sf n -2 


where 


_ 1 - a n + San-! 
1 — t + S 


(29) 


(30) 


Substituting (23) in (22) gives the final form of M n as 


M n 


f a n 


tin— 1 
—Sa n -2 


Vne N 


where 


a n = ra n -1 — Sa n -2 V* £ N 


(24) 


(25) 


with initial conditions do = 1 and a_i = 0. Apart from 
the iterative definition, it is also possible to explicitly 
determine a n in terms of known quantities. 

In terms of r and S , a n is given as 


n 
2 

a n =J2 (- 1 )™ n ~ m Cm s m T n ~ 2m Vn > 0 (26) 

m —0 


and a_i = 0. Here [•] is the greatest integer function and 
n C r is the coefficient of x r in the binomial expansion of 
(1 + x) n . Using the properties of n C r , it can be shown 
the a n as expressed in (26) satisfies (25). However as the 
proof is lengthy, it is given in Appendix C. 

We can also obtain another representation of a n if the 
results are expressed in terms of the eigenvalues of M, 


W,2 



(27) 


To do so, we put r = Ai + A 2 and S = A 1 A 2 in (21). 
We then decompose M as M = UDU -1 where D is the 
diagonal matrix with Ai and A 2 on its diagonals and U 
is the matrix with the eigenvectors of M as the columns. 
Then M n = UD n U ~ 1 , which when computed explicitly 
gives 


M n 


1 ( A” +1 -A? +1 Aj — A” \ 

A 2 - Ai \A” +1 A 2 - A” +1 A! A?A 2 - A^A J 


The detailed proof of this result is given in Appendix E. 
Note that the proof is independent of the explicit form 
of a n and uses only the propagation rule (25). 


VII. COMMENTS ON SPEED OF THE 
ALGORITHM 

Although results in Section-5 demonstrate the use of 
the algorithm described in this article, the effectiveness 
of the proposed algorithm can be gauged better in higher 
dimensional systems where computation of existence and 
stability conditions becomes computationally intensive 
due to large orders of matrices involved. In this section, 
we show the efficiency of the proposed method against 
the traditional methods as the dimensions of the system 
increases. To do this we compare the times required by 
the proposed method to compute powers of the matrix 
M l /r for various dimensions against the times required 
by other frequently used methods for the same computa¬ 
tion. 

First let us look at the computational complexity of 
the proposed algorithm. As noted erlier, the efficiency of 
the proposed method lies in the structure of the matrices 
Ml/r- We note from (17) that 

WX, = r^-o'-D = r ii(n+1) _ ((j . +1) _ 1) = [M n+ \. +1 . 

(31) 

This implies that in V-dimensions, the first N — 1 
columns of M n is identical to the last N — 1 columns 
of M n+1 . Hence, to compute M n+1 from M n , we need 
to compute only one new column of the matrix. The con¬ 
struction of this new column is described by (15). Com¬ 
putation of each term of the column requires N multipli¬ 
cations. As there are N elements in the row, the total 
number of multiplications to be performed to raise the 
V-dimensional matrix M to the power n is N 2 n. Hence, 
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the proposed algorithm has the computational complex¬ 
ity of the order 0(N 2 n). 

Two of the most common methods for computing the 
powers of matrices are brute-force matrix multiplication 
and matrix diagonalisation. Brute-force matrix multipli¬ 
cation simply multiplies the matrices successively; hence 
making computation of n th power of an iV-dimensional 
matrix an 0(N 3 n ) process. On the other hand, ma¬ 
trix diagonalisation seeks to diagonalise the matrix M 
as M = UDU~ 1 and then compute the n th power as 
M n = UD n U^ 1 . While the method seems elegant, it in¬ 
volves diagonalisation of the matrix; making it an 0(N 3 ) 
process [21]. Even with more sophisticated matrix mul¬ 
tiplications like Strassen algorithm [22] or Coppersmith- 
Winograd algorithm [23-25], one obtains a order com¬ 
plexity of 0(N q n) with q > 2. Comparing the order 
complexity of the proposed method (which is 0(N 2 n)) 
with that of brute force matrix multiplication or matrix 
diagonalisation, we see that the proposed algorithm is 
better than the other methods at least for large enough 
matrix dimension. 

However, to get an estimate for realistic dimensions 
and powers, we computed the powers of random ma¬ 
trices of the form (10) on a computer using the three 
algorithms: a) the proposed algorithm, b) matrix di¬ 
agonalisation and c) brute force matrix multiplication 
to compare their run-times. Figure 2 shows the time 
taken by the different algorithms on the same computer 
to compute the powers of random matrices of the form 
(10). As can be evidently seen, the efficiency of the pro¬ 
posed method, when compared to the other algorithms 
increases with increase in the dimension of the matrices 
when the power to which the matrices are raised is kept 
constant. When the matrix power is increased keeping 
the dimension of the matrix fixed, we see that the gain 
provided by the proposed algorithm as compared to ma¬ 
trix diagonalisation reduces as n is increased. However, 
the proposed algorithm performs better than matrix di¬ 
agonalisation for a significantly large range of n. Only 
for very high powers of matrices is matrix diagonaisation 
a better algorithm than the proposed one. 

Hence the proposed algorithm is seen to perform bet¬ 
ter than the competing algorithms for high dimensional 
systems when we are interested in sufficiently lower pow¬ 
ers of the matrix. Physically, this corresponds to pe¬ 
riodic orbits of sufficiently low periodicity in high di¬ 
mensional phase space. Non-smooth dynamical systems 
with high dimensionality occur in many real-life electri¬ 
cal, electronic and robotic systems where the dimensions 
of the phase space can well go over 50 due the presence of 
many components. The proposed algorithm might hence 
be applied to such dynamical systems to increase com¬ 
putational efficiency. 



Dimensions of the Matrix (N) 


(a) Power computed, m = 10 



(b) Dimension of the Matrix, N = 50 


FIG. 2: Comparison of the proposed algorithm of 
raising matrices to powers with the pre-existing 
methods. Plots show the time taken to perform 100 
typical calculations by each of the algorithms. For each 
dimension and power, 1000 sets of 100 random matrices 
of the form (10) each were taken and the total time 
taken for computation was averaged over the 1000 sets. 


VIII. CONCLUSION 

In this article we developed a faster and a more el¬ 
egant technique to compute the existence and stability 
conditions for periodic orbits of the form L m R n in an 
arbitrary dimensional piecewise linear continuous map. 
The technique is based on easier computation of powers 
of N x N matrices in their normal form. Due to the 
structure of the matrices involved, it was found that the 
elements of the resulting matrix were interrelated; and 
in order to compute the n th power of the matrix, only 
N out of the N 2 elements need to be computed. These 
N elements in turn can be obtained as simple sequences 
defined iteratively. Once the powers of the matrices are 
computed, they can be substituted in the generic expres¬ 
sions of existence and stability of orbits to obtain the 
regions in parameter space where they exist. We also 
apply the technique developed to 3 dimensional systems 
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and obtain the regions where L n R orbits exist. 

Moreover, in the special case of 2 dimensional matri¬ 
ces, further simplifications were made. Notably, explicit 
expressions for the terms of the sequence in terms of the 
given parameters and eigenvalues of the matrices were ob¬ 
tained. This allows for a direct evaluation of any power 
of the 2x2 normal form matrix without computing the 
intermediate powers. 
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Appendix: Derivations and Motivations 


Appendix A: Finding M n for a Matrix M in the Normal Form 


Theorem A.l Let M a matrix of the form (10). Then the elements of M n are given as 

[M n ] • • = D,„_(*_!) Vn e N (Al) 

where [A] i . is the element of A corresponding to i th row and j th column and F ,j is defined in (13) and (15). 

Proof Let us assume the most general form of M n 

(oft ••• 0§\ 


for some n > 1. Then 


M n = 


a( n ) 
\"N, 1 


,(n) 


,n/ 


M n+1 = 


/n(n) n(n) 
t7 l,l "1,2 
Jn) Jn) 
" 2,1 u 2,2 


y N,N 


g(n) \ 

°1,N 

fl(") 
°2 ,N 


Pi 

~P2 

\N -1 


/^n+i) 

/}( n +l) 

^1,2 

■ ot n i] \ 

^2,1 

/i(n+l) 

^2,2 

^(n+l) 
u 2 ,N 

ya(n+l) 

\^AT,1 

/^(n+1) 

U N, 2 

fi( n + 1) ; 
U N,N / 


2,1 


(") \ 
l.JV-l 

q( n ) 

y 2,N-l 


°\ 

0 


VA <1 ■■■ Vi-D"-' PN «■■■«) 

/(eLi(-i••• » 


(A2) 


Hence 


and 


5( ” +1) = -^<i<N,n>l 


iV 


,(n) . 


/e—1 


^” +1) = ^S -1 : 1 < i < N, l<j<N,n>l. 

We now use (A4) iteratively to obtain 


Now if we define 


then 


6 ft +1) = = • • • = °“ 2)) : 1 < * < N, 1 < j < N, n > 1. 


dlT 1] = D,n+i : 1 < i < N 


0g +1) = 4”r (i_2)) = D,n- 0 - 2 ) : 1 < * < IV, 1 < j < N, n > 1. 


Combining (A6) and (A7); and replacing n + 1 by n we get 


0$ = r i]n _ (j -_i) : 1 < i,j < N,n > 2. 


(A3) 

(A4) 

(A5) 

(A6) 

(A7) 

(AS) 


Now, by definition 


N 


r i ,„ +1 = < n ! +1) = ]T(-i) fe -V0. 


(") 

i,k 


k =1 


: 1 < i < iV,n > 1. 


(A9) 
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Using (A8) in the right hand side gives us 

N 

r i>n+1 = fl£ +1) = ^(—i) fc - 1 p fc r ij „_ (fc _ 1) : 1 < * < N,n > 1. 

fc=l 

Finally replacing n + 1 by j, we have an iterative relation for Tij as 

N 

T u = 5I(-l) fc "Vfc: 1 < i < JV, j > 2. 

k=l 


Now note that substituting n = 1 in (A8) gives 


9?) =T i>2 -j :1 <i,j<N. 


However by dehnition of , we have 


9^ = [M 1 ]. J =M i , j : 1 < i,j < -V. 


Hence, 


r il2 -i=M w • 1 < i,j < N 


which on replacing 2 — j by j yields 


Uj = Mi, 2 -j ■■ 1 < i < N, -N + 2 < j < 1 


(A10) 

(All) 

(A12) 

(A13) 

(A14) 

(A15) 


Appendix B: Motivating the Form of a n 

In this Appendix, we give the motivation for obtaining the fundamental sequence in the form 


ctn=J2 n ~ mc m <5 m T n - 2m Vn > 0 

m —0 

and a_i = 0. In other terms, we try to understand, how the matrix 

M = 

yields 

M n = 


r 1 

—8 0 


Un O’n—l 

8&n — l 8(l n —2 


Vn £ N 


with a n dehned in (Bl). 

For an matrix general 2x2 matrix M, if we assume 


_ I ®n b n 

. Cn d n 


then the sequence 


or 


ai 6i 
ci di 


M, M 2 ,M 3 ,... ,M r 


a-i fa 3 b 3 \ fa n b r 

C 2 d 2 ) ’ lc 3 d 3 ) \c n d r 


(Bl) 

(B2) 

(B3) 

(B4) 

(B5) 

(B6) 
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is a set of four sequences: {a n } , {b n } , {c n } and {d n }. 

The aim of the section is to show that these four sequences are restricted by constraints that allow the matrix to 
be expressed in terms of a single sequence. 

To show this, we assume that (B4) holds for the matrix defined in (B2). Then 

M n+1 = M n .M 


fdn- f-1 ^n+l\ _ fd n b n \ { T l\ 

yC n _)_i d n -\-l J yC n d n J y 6 0 J 

{ 1 ^n+l\ _ fTd n db n UnA 

yCn+1 d n -\-\ J \TC n Sd n C n J 


Hence 


or 


bn+l d n 

dn -\-1 — Cn 


bn — dn— 1 


(B7) 


and 


dn — C n —\ 


d n +i — Ta n db n 


(B8) 


(B9) 


C n+ 1 = TC n - Sd n . (Bio) 

Substituting (B7) in (B9) and (B8) in (BIO), we get 

a n + i = ra n - Sa n -1 (Bll) 


Cn+l — TCn dCn-1- 


Hence, at the current stage 


M n 


f CL n &n—l 
1 


(B12) 


(B13) 


with a n and c n satisfying (Bll) and (B12) respectively. 
Moreover, as we know 


therefore 

ai = r, ao = 1, ci = —(5, cq = 0. (B14) 

Using (Bll) and (B12) in conjugation with the initial conditions in (B14), can write the complete sequence of a n 
and c n . The first few terms are shown below. 
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0,1= T 

a 2 = T 2 - 6 

03 = t 3 — 2t5 

ai=T A — 3t 2 S + S 2 

05 = r 5 — 4 t 3 (5 + 3 t5 2 

ae = r 6 — 5 t 4 <5 + 6t 2 <5 2 — S 3 

ar = t‘ — 6r S (5 + 10r 3 (5 2 — 4r<5 3 

o 8 = t 8 - 7t 6 S + 15t 4 S 2 - 10t 2 6 3 + 6 4 

a g = T 9 - 8t 7 6 + 21 t 5 S 2 - 20 t 3 S 3 + 5 t6 4 
It may be noted that if written in the appropriate form. 

c n — do, t 


ci = —S 
ci = —<5[r] 
c 3 = — 6 [t 2 - <5] 

C4 = — <5[r 3 — 2tS\ 

C 5 = — S[t 4 — 3t 2 5 + 5 2 \ 
cq = — <5[r° — 4 t 3 S + 3t6 2 } 

C 7 = —<5[r 6 — 5 t 4 S + 6t 2 S 2 — <5 3 ] 
eg = — 5[t‘ — 6t 5 S + 10 t 3 8 2 — 4r(i 3 ] 

eg = — <5[r 8 — 7 t 6 S + 15r 4 ^ 2 — 10r 2 <5 3 + d 4 ] 
it becomes clear that 

Vn e N. (B15) 


In order to extend the result to eg, we define a_i = 0 which allows us to write M n as 


M n 


( a n 
y 3o n —i 


Q"n —1 

—Sa n -2 


(B16) 


where 


a n +i = Ta n — Sa n - 1 Vn > 0 (B17) 

and a_i = 0. 

To obtain the analytic expression for a n , we note the following in expansions given earlier, when terms under 
summation for a particular n value are arranged in the ascending order of powers of 5, 

• There are [^-] + 1 terms in the series. 

• Powers of S start from 0 and increase in steps of 1. 

• Powers of r start from n and decrease in steps of 2. 

• Sign of each term alternates between plus and minus starting from a plus 

• A numerical coefficient precedes each term. 

Hence a n can be written as 

[§] 

a n =J2 (-!) m Vm,n 6 m T n - 2m (B18) 

m —0 

A look at the 7] m ^ n values in Table I reveal that 

Vm,n = Vm,n — 1 4” Vm — l,n—2 (B19) 

which seems similar to the properties of the binomial coefficients. In fact, 

Vm,n = n ~ m C m (B20) 

and the structure of the Pascal’s triangle might be evidently seen in the table. 


Appendix C: The Progression Rule of Fundamental Sequence 


Lemma C.l The n th term of the sequence defined in (26) is related to its two preceding terms by the relation 

On = ra n —i do n - 2 - 


(Cl) 
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TABLE I: List of all rjm^n values 



0 12 3 4 

l 

1 

2 

1 1 

3 

1 2 

4 

1 3 1 

5 

1 4 3 

6 

15 6 1 

7 

1 6 10 4 

8 

1 7 15 10 1 

9 

1 8 21 20 5 


Proof For n = 0, 

a n +i = cl\ = t 
O n — a 0 — 1 
&n— 1 = 0—1 = 0 

Hence (Cl) is true for n = 0. 

For the other n, we prove it separately for even and odd n. 

If n is even, then 


n 

n 

n — 1 

1 -2-1 

~ 2’ 

2 


= - - 1, 
2 


n + 1 


n 

2 ' 


Now, 


n — 1 


_ e _ \ A ({ i\m rm n-2m\ e \ A // i\m n—1 —rm n-l-2m\ 

ra n - da n _i = r 2^ ((-1) C m d r J - d 2_^ (( _1 ) 6 m d r J 

ra—0 m=0 

n n 

2 2 _1 

= 53 ((” 1 ) m n ~ mC m 5 m T n+1 ~ 2m ) ^ 53 ((“I)™ n ~ X ~ m C m S^t™- 1 - 2 ™) 


m =0 


m —0 


n 

2 


n 

2" 1 


= T n+1 + 53 ((-l) m n ~ m c m 6 m T n+1 ~ 2m ) - Y ((-!) m n ~ l ~ m C m J m + 1 r n " 1 ' 2m ) 

m =1 m—0 

n n 

2 2 

= T n+1 + 53 ((- 1 ) m n ~ m Cm S m T n+1 - 2m ) + Y ((- 1 )” 1 n ~ m Cm -1 (TV"* 1-2 "*) 


m=l 

n 

2 


m=l 


= r n+1 + 53 ((-i) m ( n ' m c m + n ~ m c m _ i) r T n+1 - 2m ) 

ra=l 

n 

2 

= T n+1 + 53 (( -1 ) m n+1-m C' m ^ m T- n+1_2rn ) 


m =1 


n + 1 


= 53 (( _1 ) m ” +1_m Cm ,5 m T- rl + 1_2rn ) 


m =0 


(C2) 


— On+l- 
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If n is odd then 


Now, 


w 

n — 1 

n — 1 

n — 1 

n + 1 

71 + 1 

L 2 J 

2 ’ 

2 

2 ’ 

2 

2 


(C3) 


~ \ (( i \m n—m/~i cm n-2m\ c \ A // i\m n —1—em_n-l-2m\ 

ra n - da n _i = r 2^ ((-1) C m o r ) - d ^ ((-1) C m o r ) 


m —0 

a — 1 


m —0 


= ^ ((-l) m n ~ m c m s m T n+1 ~ 2m ) - Y ((- 1 ) 




') 


m—0 


m —0 


n — 1 „ n -\-1 


= T n+1 + ^ ((-I)"" ,5 m T "+ 1 - 2m ) - Y ((-I)"" ,5rn+l r n-l-2m^ _ (_!)+*-£ 

m= 1 m=0 

71 — 1 71-1 

2 2 

= T n+1 + Y ((“ 1 ) m n ~ m Cm 5 m T n+1 ~ 2m ) + Y ((- 1 ) 1 " n ~ m Cm- 1 <TV” +1 - 2ro ) + (-1)^5^ 

m=l m— 1 

71-1 

2 

= T n+1 + Y ((-!) m ( n ~ m Cm + n ~ m C m _ 1 ) ( 5 m T "+ 1 - 2m ) + (-l)^^ 

m=l 

71-1 

2 

= T n+1 + ^ ((-I)™ n+1 ~ m C m S m T n+1 - 2m ) + (-1 


m =1 


ri + ll 


= Y ((-!) m n+1 ~ m Cm S m T n+1 ~ 2m ) 

m—0 
= CLn -\-1 • 

Hence the result is true for all n£ff. 


Appendix D: M n in Terms of Eigenvalues 


Theorem D.l Let 


with the eigenvalues 


then 



Al,2 



(Dl) 


(D2) 


M n 


Q"ri — 1 

Sa n —i Sa n — 2, 


Vn G N 


(D3) 


where 


a n 


A" +1 


A£ +1 


Ai — A2 


(D4) 


Proof Let Ai and A 2 be the eigenvectors of M. Substituting the trace r = Ai + A 2 and determinant 8 = AiA 2 , we 
have 


M = 


Ai + A 2 1 
—AiA 2 0 


(D5) 
















14 


Simple calculation of shows that the eigenvector corresponding to Ai is and that corresponding to A 2 is 

\-A 2 / 


. —Ai 


|. Hence, we may construct a matrix U with eigenvectors as columns, 


and hence 


Substituting the values, we get 


M n = 


U = 


1 1 

. —A 2 —Ai 


and a diagonal matrix D with the eigenvalues 


such that 


D = 


Ai 0 
, 0 A 2 


M = UDU- 1 . 


M n = UD n U~ 1 . 


1 1 Ai 0 W 1 1 

, — A 2 — Ai/ \ 0 A 2 / A 2 — Ai 


(D6) 

(D7) 

(D8) 

(D9) 


1 f 1 1 ) 0 ] f-Ai -1 

A 2 — Ai l — A 2 — Ai J l 0 A 2 / l A 2 1 


1 1 \ I —A” +1 -A” 


A 2 Ai \ —A 2 —Ai I \ A 


n+1 


which can be recast as 


with 


A«+i _ \ n + 1 


A£-A? 


A 2 -Ai l A^ +1 A 2 -A” +1 Ai A”A 2 — A 2 Ai 


M n = 


O'n tin— 1 

ScL n —1 S(l n —2 , 


— 


^n+1 _ ^n+1 

Ai — A 2 


(DIO) 


(Dll) 


Appendix E: The Form of 1 
Corollary E.l The sum of the geometric progression 


is given as 


= I + M + M 2 + ... + M n 


fn fn—1 


(El) 
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where 


1 d n “b Sd n — 1 

1 — T + 5 


Proof Using the formula for sum of GP, we can write </>„ as 




I — M n 
I -M ' 


Hence 




O'n— 1 

1 + Sa n -2 


(E2) 


(E3) 


1 — r + <5 \ — S 1 — 



1 ftn (in-1 

5a n -i 1 + Sa n —2 , 


1 

1 — r + (5 


1 ei n Sd n —i 1 o, n —i ~b Sd n —2 

.-5 (1 - an) + Sdn- 1 (1 - r) Sdn-! + (1 - t) (1 + Sd n - 2 ) 


1—a n +<5a n i 


1— a n —i +<5a n 2 


1—r+<5 1—r+£ 

—(5(1—Q n )+<5a rt -i(l—r) <5an-i + (l —r)(l+^a n 
1—r+<5 1 —t+<5 


Currently, <f> n is of the form 


with 


and 


(7l = 


ct 2 = 


</>„=(/" fn ~ 1 ) 

\(7i cr 2 y 

-(5 (1 - a n ) + (5a n _i (1 - r) 
1 — r + (5 

fe„_i + (1 - t) (1 + Sd n - 2 ) 

1 — T + S 


Using lemma C.l, we simplify o\ and cr 2 as 


(7l = 


-S (1 - a n ) + Sa n - 1 (1 - t) 

1 — t + S 

— S + Sd n “b Sd n — i — STd n —i 
1 — t + S 

-S (1 - d n - d n -l + Td n - 1) 

1 — T + 5 

— S (1 — d n — On~l + dn + Sd n — 2 ) 
1 — T + S 

-S (1 - a n -1 + Sa n - 2 ) 

1 — T + S 
= Sfn—1 • 


(E4) 


(E5) 


(E6) 
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And 


Therefore, 


_ fan -1 + (1 - t)(1 + Sa n - 2 ) 

1 — t + S 

_ fan -1 + 1 + Sa n - 2 T T Sa n —2 

1 — T + S 

_ 1 - r + (5 (a ra _i + a n -2 - 'ra rl _ 2 ) 

1 — T + (5 

1 — r + (5 (a„_i — fa n - 3 ) 

1 — r + S 

_ 1 ~ T + (5 (-1 + 1 + Ort-i - fart- 3 ) 
1 — r + (5 

_ 1 - r + d + 8 (-1 + a n -i - Sa n - 3 ) 
1 — r T d 

_ 1 — t + (5 — (5 (1 — a n _i + fan- 3 ) 

1 — r + (5 

_ 1 (5(1 — a ra _i + fan_ 3 ) 

1 — r + d 

= 1 - (5/„_ 2 . 


And 


Therefore, 



(E7) 










